home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / gs_iter.pro < prev    next >
Text File  |  1997-07-08  |  7KB  |  196 lines
  1. ;$Id: gs_iter.pro,v 1.7 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       GS_ITER
  8. ;
  9. ; PURPOSE:
  10. ;       This function solves an n by n linear system of equations
  11. ;       using Gauss-Seidel iteration.
  12. ;
  13. ; CATEGORY:
  14. ;       Linear Algebra.
  15. ;
  16. ; CALLING SEQUENCE:
  17. ;       Result = GS_ITER(A, B)
  18. ;
  19. ; INPUTS:
  20. ;       A:      An N by N array of type: int, float, or double.
  21. ;
  22. ;       B:      An N-element vector of type: int, float, or double.
  23. ;
  24. ; KEYWORD PARAMETERS:
  25. ;       CHECK:    An integer value of 0 or 1 that denies or requests
  26. ;                 checking A for a diagonally dominant form.
  27. ;                 CHECK = 0 (the default) results in no checking.
  28. ;                 CHECK = 1  Checks A and reports if it does not
  29. ;                            meet the required condition. This is
  30. ;                            just a warning. The algorithm will
  31. ;                            proceed on the chance it may converge.
  32. ;
  33. ;       LAMBDA:   A scalar value in the range: [0.0, 2.0]
  34. ;                 This value determines the amount of 'RELAXATION'.
  35. ;                 Relaxation is a weighting technique that is used
  36. ;                 to enhance convergence.
  37. ;                 1) LAMBDA = 1.0 (the default) no weighting.
  38. ;                 2) A value in the range  0.0 <= LAMBDA < 1.0  improves
  39. ;                    convergence in oscillatory and non-convergent systems.
  40. ;                 3) A value in the range  1.0 < LAMBDA <= 2.0  improves
  41. ;                    convergence in systems known to converge.
  42. ;
  43. ;       MAX_ITER: The maximum number of iterations allowable for the
  44. ;                 algorithm to converge to the solution. The default 
  45. ;                 is 30 iterations.
  46. ;         
  47. ;       X_0:      An N-element vector that provides the algorithm's 
  48. ;                 starting point. The default is [1.0, 1.0, ... , 1.0].       
  49. ;
  50. ;       TOL:      The relative error tolerance between current and past
  51. ;                 iterates calculated as:  ABS( (current-past)/current )
  52. ;                 The default is 1.0e-4.
  53. ;
  54. ; SIDE EFFECTS:
  55. ;       Upon output A and B are divided by the diagonal elements of A.
  56. ;       Integer inputs are converted to floats.
  57. ;       Note: These SIDE EFFECTS have been removed for IDL v5.0.
  58. ;
  59. ; RESTRICTIONS:
  60. ;       The equations must be entered in a DIAGONALLY DOMINANT form
  61. ;       to guarantee convergence.
  62. ;       A system is diagonally dominant if it satisfies the condition:
  63. ;                   abs(A(row,row)) > SUM(abs(A(row,col)))
  64. ;       where SUM runs col=1,N excluding col = row and A is in row major.
  65. ;       This restriction on A is known as a sufficient condition. That is,
  66. ;       it is sometimes possible for the algorithm to converge without the
  67. ;       condition being satisfied. But, convergence is guaranteed if the
  68. ;       condition is satisfied.
  69. ;
  70. ; EXAMPLE:
  71. ;       Define an array (a) in a non-diagonally dominant form.
  72. ;         a = [[ 1.0,  7.0, -4.0], $
  73. ;              [ 4.0, -4.0,  9.0], $
  74. ;              [12.0, -1.0,  3.0]]
  75. ;       And right-side vector (b).
  76. ;         b = [12.0, 2.0, -9.0]
  77. ;       Compute the solution of the system, ax = b.
  78. ;         result = gs_iter(a, b)
  79. ;       Note: This example fails to converge, because A is not in
  80. ;             diagonally dominant form.
  81. ;
  82. ;       Reorder the array given above into diagonally dominant form.
  83. ;         a = [[12.0, -1.0,  3.0], $
  84. ;              [ 1.0,  7.0, -4.0], $
  85. ;              [ 4.0, -4.0,  9.0]]
  86. ;       Make corresponding changes in the ordering of b.
  87. ;         b = [-9.0, 12.0, 2.0]
  88. ;       Compute the solution of the system, ax = b.
  89. ;         result = gs_iter(a, b)
  90. ;
  91. ; PROCEDURE:
  92. ;       GS_ITER.PRO implements the Gauss-Seidel iterative method with
  93. ;       over- and under- relaxation to enhance convergence.
  94. ;     
  95. ; REFERENCE:
  96. ;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
  97. ;       Erwin Kreyszig
  98. ;       ISBN 0-471-55380-8
  99. ;
  100. ; MODIFICATION HISTORY:
  101. ;       Written by:  GGS, RSI, April 1993
  102. ;       Modified:    GGS, RSI, February 1994
  103. ;                    1) Format keyword is no longer supported. The matrix
  104. ;                       should be supplied in a row major format. 
  105. ;                    2) The input/output parameter X has been removed. The 
  106. ;                       algorithm's initial starting point is an n-element
  107. ;                       vector of 1s. The keyword X_0 has been added to 
  108. ;                       override the default.
  109. ;                    3) GS_ITER is now called as a function, x = gs_iter( ). 
  110. ;       Modified:    GGS, RSI, April 1996
  111. ;                    The input arguments are no longer overwritten.
  112. ;                    Added DOUBLE keyword. Modified keyword checking and use 
  113. ;                    of double precision.
  114. ;-
  115.  
  116. FUNCTION GS_ITER, A, B, Check = Check, Lambda = Lambda, Max_Iter = Max_Iter, $
  117.                         X_0 = X_0, Tol = Tol, Double = Double
  118.  
  119.   ON_ERROR, 2  ;Return to caller if error occurs.
  120.  
  121.   TypeA = SIZE(A)
  122.   TypeB = SIZE(B)
  123.  
  124.   if TypeA[TypeA[0]+1] lt 2 or TypeA[TypeA[0]+1] gt 5 then $
  125.     MESSAGE, "Input array (A) must be integer, float, or double."
  126.  
  127.   if TypeB[TypeB[0]+1] lt 2 or TypeB[TypeB[0]+1] gt 5 then $
  128.     MESSAGE, "Input vector (B) must be integer, float, or double."
  129.  
  130.   if TypeA[TypeA[0]-1] ne TypeA[TypeA[0]] then $
  131.     MESSAGE, "Input array must be square."
  132.  
  133.   ;If the DOUBLE keyword is not set then the internal precision and
  134.   ;result are determined by the type of input.
  135.   if N_ELEMENTS(Double) eq 0 then $
  136.     Double = (TypeA[TypeA[0]+1] eq 5 or TypeB[TypeB[0]+1] eq 5)
  137.  
  138.   ;Set default values for keyword parameters
  139.   if KEYWORD_SET(Lambda)   eq 0 then Lambda = 1.0
  140.   if KEYWORD_SET(Max_Iter) eq 0 then Max_Iter = 30
  141.   if KEYWORD_SET(X_0)      eq 0 then X_0 = REPLICATE(1.0, TypeA[1]) 
  142.   if KEYWORD_SET(Tol)      eq 0 then Tol = 1.0e-4
  143.  
  144.   ;Diagonal elements of input matrix.
  145.   Diag = A[LINDGEN(TypeA[1]) * (TypeA[1]+1)]
  146.  
  147.   if KEYWORD_SET(Check) ne 0 then begin
  148.     Sum = TOTAL(ABS(A), 1, Double = Double) - ABS(diag)
  149.     caution = WHERE(Sum ge ABS(Diag), Count)
  150.     if Count ne 0 then begin
  151.       PRINT, "Input matrix is not in Diagonally Dominant form." & $
  152.       PRINT, "Algorithm may not converge."
  153.     endif
  154.   endif
  155.  
  156.   ;Precondition inputs.
  157.   ;Divide the rows of A and B by the diagonal elements of A.
  158.   if Double eq 0 then begin
  159.     AA = A / (REPLICATE(1.0, TypeA[1]) # Diag)
  160.     BB = B / (Diag + 0.0)
  161.     X_0 = FLOAT(X_0)
  162.   endif else begin
  163.     AA = A / (REPLICATE(1.0d, TypeA[1]) # Diag)
  164.     BB = B / (Diag + 0d)
  165.     X_0 = DOUBLE(X_0)
  166.   endelse
  167.  
  168.   Cond = 0
  169.   Iter = 0
  170.  
  171.   ;Begin the computational loop and continue WHILE
  172.   ;the number of iterations is less than max_iter
  173.   ;AND the relative error between iterations is
  174.   ;greater than tol.
  175.   while(Iter lt Max_Iter and Cond eq 0) do begin
  176.     Cond = 1
  177.     Iter = Iter + 1
  178.   ;Formulate x_0 as the row vectors of A.
  179.   for k = 0, TypeA[1]-1 do begin
  180.     xLast = X_0[k]
  181.     X_0[k] = Lambda * (BB[k] - (TOTAL(X_0*AA[*,k],1, Double = Double)) + $
  182.              (AA[k,k] * X_0[k])) + (1.0 - lambda) * xLast
  183.     if Cond eq 1 and X_0[k] ne 0.0 then begin
  184.       Error = ABS((X_0[k] - xLast) / X_0[k])
  185.       if Error gt Tol then Cond = 0
  186.     endif
  187.   endfor
  188.   endwhile
  189.  
  190.   if Iter ge Max_Iter and Cond eq 0 then $
  191.     MESSAGE, "Algorithm failed to converge within given parameters."
  192.  
  193.   RETURN, X_0
  194.  
  195. END
  196.